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Abstract 



Given a sequence of N positive real numbers {ai, da, . . . , cln}, the number par- 
titioning problem consists of partitioning them into two sets such that the absolute 
value of the difference of the sums of cij over the two sets is minimized. In the case 
that the Oj's are statistically independent random variables uniformly distributed 
in the unit interval, this NP-complete problem is equivalent to the problem of find- 
ing the ground state of an infinite-range, random anti-ferromagnetic Ising model. 
We employ the annealed approximation to derive analytical lower bounds to the 
average value of the difference for the best constrained and unconstrained parti- 
tions in the large N limit. Furthermore, we calculate analytically the fraction of 
metastable states, i.e. states that are stable against all single spin flips, and found 
that it vanishes like N~ 3 / 2 . 



Short Title: number partitioning problem 
Physics Abstracts: 87.10.+e - 64.60.Cn 



1 



1 Introduction 



The importance of the study of complex optimization problems which involve 
quenched, random, frustrated functions of many variables, as well as the major 
role that statistical mechanics can play in that study, have been pointed out by 
Anderson more than ten years ago \%. Since then, standard statistical mechanics 
techniques have been applied to the probabilistic analysis of several classical com- 
binatorial optimization problems, such as the graph partitioning problem Q, the 
traveling salesman problem |j, || , the knapsack problem JB], [| fjj , and the satisfiabil- 
ity problem ^ ||, [l(J , to mention only a few. In fact the well-established statistical 
mechanics methods to characterize ground states (global minima) and metastable 
states (local minima) of spin glass models can be readily adapted to the study of 
optimization problems [Tlfl . 

In this paper we study the number partition problem (NPP) which is stated 
as follows. Given a sequence of real numbers {ai, ai 1 ■ ■ ■ , cln}, the NPP consists of 
partitioning them into two disjoint sets A\ and Ai such that the difference 



E 

Oj£.4i a,G^2 



is minimized. Alternatively, we can search for the Ising spin configurations s = 
(si, . . . , Sjv) that minimize the energy or cost function 

N 

E(s)= |5>i*il. (2) 
j=i 

where Sj — 1 if a.j G A\ and Sj = — 1 if ctj € A^- We can consider also the problem 
of constrained partitions, in which the difference between the cardinalities of sets 
A\ and Ai is fixed, i.e., 

1 N 

m=- ^Sj. (3) 

3=1 

The NPP may be viewed as the practical problem of finding the fairest way to 
partition a set of N objects j = 1, 2, . . . , N, each of which of value aj, between two 
persons. Despite its simplicity, the NPP was shown to belong to the NP-completc 
class, which basically means that there is no known deterministic algorithm guaran- 
teed to solve all instances of this problem within a polynomial time bound [|l2| . The 
fact that the NPP is frustrated can easily be understood by squaring equation (||) , 
so that the problem of minimizing E becomes then the one of finding the ground 
state of the infinite-range, random anti-ferromagnetic Ising Hamiltonian [fill 



i j>i 

Thus we note that the problem of finding the ground state of (0) is NP-complete. 
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Although zero-cost solutions of the NPP may be of some value to cryptography 
[ H , the interest in this problem stems mainly from the remarkable failure of the 
stochastic heuristic simulated annealing jl^, [l6j to find good solutions to it, as com- 
pared with the solutions found by deterministic heuristics jlTj]. In fact, the reason 
for that failure is that the usual strategy of exploring the space of configurations 
{s} through single spin flips leads to changes of energy that are typically of order 
1/iV, while a theoretical analysis indicates that the global minimum energy is of 
order t/N 2~ n for unconstrained partitions fTsfl . It is interesting to note that a 
very simple deterministic heuristic, the differencing method of Karmakar and Karp 
[ [l9[ , can find with high probability solutions whose energies are of order 1 /N a log N 
for some a > 0. More recently, it has been shown that the performance of simu- 
lated annealing can be greatly improved and even surpass that of the differencing 
method by employing different representations for the problem |pof . 

In this work we employ the annealed approximation [pl| , |22f to derive rigorous 
lower bounds to the average value of the difference or energy for the best constrained 
and unconstrained partitions. For constrained partitions, we show that the average 
optimal energy is extensive for m > \/2 — 1 and we calculate it exactly in this 
regime using the self-averaging property of the free energy density. The theoretical 
predictions are compared with numerical estimates for the optimal energy obtained 
through the exhaustive search of the configuration space for N < 24. Furthermore, 
we calculate analytically the average number of minima in the 1-swap neighborhood 
and estimate their typical energy. A minimum in the 1-swap neighborhood is a state 
that has lower energy than all the N states that differ from it by a single spin only 
[ p"7[ . Similarly to previous studies of the NPP jlT], [l8| |2(| , we will consider the case 
where the Oj's are statistically independent random variables uniformly distributed 
in the unit interval. 

The remainder of this paper is organized as follows. In section 2 we describe 
the annealed approximation and calculate the lower bounds to the average value of 
the optimal energy. In section 3 we present the calculation of the average number of 
local minima in the 1-swap neighborhood. In section 4 we discuss our main results 
and present some concluding remarks. In particular, we compare our approach 
with other theoretical studies of the NPP jl3| . In the appendix we present the 
details of the self-averaging calculation of the average optimal energy in the regime 
where this quantity is extensive. 

2 Annealed approximation 

In the canonical ensemble formalism of the statistical mechanics the average value 
of the optimal energy for constrained partitions is given by 

E m = lim F m (T) — - lim T (lnZ m ) , (5) 

T~>0 T^O 
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where F m (T) is the average free energy, and Z m (T) is the partition function 

E( s y 



Z m (T) = J28\ Nm, s o ] ex P 

M 



T 



(6) 



with m = —1, —1 + 2/N, . . . , 1 — 2/N, 1. Here the summation is over the 2 N states s, 
5(k,l) is the Kronecker delta and T is the temperature. The notation (. . .) stands 
for the average over the random variables a^. The limit T — * in equation (^|) 
ensures that only the states that minimize E (s) will contribute to Z m . 

Since the average entropy S m (T) — ~dF m jdT of a system of Ising spins is 
positive at all temperatures, F m must be a decreasing function of T, so that E m — 
F(0) > F(T) for all T. Defining the annealed free energy by 

F«(T) = -T In (Z m (T)), (7) 

and using Jensen's inequality [ p3| , ln(Z TO ) > (In Z m ) , yield the following inequalities 

F r a n (T) < F m (T) < E m . (8) 

Thus, the annealed free energy calculated at any T provides a rigorous lower bound 
to E m H). Clearly, the tightest bound is given by E^ = F^(T^) where T* n is 
the temperature that maximizes F!^(T), i.e. 



dF a 
dT 



0. 



(9) 



This procedure is very useful because, in general, the annealed free energy is much 
easier to evaluate than the quenched one. 

We now proceed with the explicit evaluation of the annealed free energy. Using 
the integral representations of the Dirac and Kronecker delta functions we write 



(Z m (T)) = 



OO /'CO 



oo «/ — OO 



n 



da 



dxdx 
1 

, E 

,-=±i 



-\x\/T 



71 dm 



iNmm 



2ir 

exp [— isj (djX + m) 



(10) 



The integrals over x and ay, as well as the summation over Sj, can easily be per- 
formed yielding 

(Z m (T)) = 



dx 



2T 



-oo 2tt i + (TxY 



sin (i/2) 
x/2 



N 



71 dm 
- 2tt' 



iNrnm 



2 im-\-ix/2 _|_ ^ — im — ix/2 



N 



(ii) 
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Using the binomial theorem, the integral over rh can be readily carried out. The 
final result is simply 



ri 



where 

n = N^, (13) 

G m (y)=imy + ]n(^-Y (14) 

and we have made the change of variable y = x/2. In the limit of large N, the 
integral over y can be evaluated using the saddle-point method J24|]. Since | m |< 1, 
the saddle-point is the imaginary y s = where £ is the real solution of the equation 

m-cothC+^=0. (15) 

Hence, the function G m (y s ) — G m , where 

G m = -mC + ln^^, (16) 

is real. Finally, using Stirling's formula for the binomial coefficient we rewrite 
equation (|lj) in the limit of large ./V as 



where 



and 



"■"'l/ d-^ic;! ^^ (17 > 

1 + TO. 1 + m 1 — m.l — m ,,„. 
9m = G m — In — — In — — (18) 

G: = -l + m 2 + ^. (19) 

At this stage we can readily calculate the temperature that maximizes the 
annealed free energy. In fact, equation @ is written as 

1 4- (IT* cS 2 

In (Z m (T r ;)) + +l mU 2 = 0. (20) 
1 - (2T*C) 2 

We consider first the regime where (Z m (T m )) is of order 1. In this case, equa- 
tion ( |l7| ) implies that T m is vanishingly small, so that equation ( pp|) reduces to 
(Z m (T^J) = e _1 . Inserting this result into equation (Q) yields = T^. Hence, 



E m = ^^{l-m*)\G' m \ e- 1 "^ (21) 
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which is consistent with the assumption that is small for large N, provided that 
g m > 0. Since g rn decreases monotonically with m, from go = In 2 to g\ = — oo, 
this assumption breaks down for | m |> 0.560 where g m is negative. Henceforth we 
will assume that m > 0. 

It is instructive to consider in detail the case of even partitions (m = 0). In 
this case we find ( = 0, go = ln2, and G = —1/3 so that 

{Z Q (T))=2 N T^. (22) 

and 

EX = 2- N ^y= « 0.167 2~ N N. (23) 
4eV3 

In figures 1(a) and 1(6) we present the results of numerical experiments to estimate 
the energy of the global minima for even partitions through the exhaustive search 
in the configuration space for N < 24. In all experiments discussed in this work, the 
symbols represent the averages over 10 4 realizations of the set {a,j}. The error bars 
are calculated by measuring the standard deviation of the average optimal energies 
obtained in 25 experiments, each one involving the average of 400 realizations of the 
set {cij}. In these experiments we focus on the N dependence of the average optimal 
energy E m = (E m ), and of the ratio r m = y/a^/E m where a 2 m = (E^) - (E m ) 2 is 
the variance of the random variable E m . In figure 1(a) we show Eq as a function of 
N. The straight line shown in this figure yields the fitting E = 0.80 2~ N N. Hence, 
although the annealed bound E§ gives the correct scaling with N, it is about five 
times smaller than our numerical estimate for Eq. In figure 1(b) we show the ratio 
ro as a function of N. Interestingly, this ratio tends to 1 for large N indicating 
then that the optimal energy Eq is not self-averaging. 

In the regime where (Z m (T^)) is of order e N we find In (Z m (T^)) rj Ng m and 

w 1/2C so that 

E& = -N^. (24) 

Of course, this solution is valid only for m > 0.560 where g m is negative. We note 
that ( pi]) gives a very poor lower bond to E rn . In particular, for m = 1 we have 
Ei = N/2 while the annealed bound yields Ef = 0. Fortunately, in the regime of 
extensive E m we can use the self-averaging property of the free energy density to 
calculate E m exactly for large N (see Appendix). The final result is simply 



F - N 
rn ~ 2 



(1 + m) 2 _ 



(25) 



which is valid for m > \/2 — 1 w 0.414. Thus the annealed lower bound is also 
very poor in the region 0.414 < m < 0.560 since in this region E^ decreases 
exponentially with N, while E m actually increases linearly with TV. 
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To better appreciate the qualitative differences between the regimes of distinct 
scalings with N, we present in figure 2 the numerical estimates for E m as a function 
of m for N — 24. The existence of two different regimes of scaling with TV, as well 
as the very good agreement with the theoretical predictions for m > 0.414, are 
apparent in this figure. A noteworthy feature of our numerical estimate for E m 



shown in the inset is that, in contrast to the annealed lower bound (21), the even 
partitions (m = 0) do not give the lowest energy. We have verified that this result 
holds for smaller values of N as well. Furthermore, there seems to occur a rather 
abrupt transition at m « 0.25 as indicated by the large error bar and for the change 
of almost three orders of magnitude in E rn . Although it would be very interesting 
to study these results more carefully for larger N, we are not aware of any efficient 
heuristic to solve the NPP for constrained partitions. In particular, we note that 
the differencing method |19 applies only to unconstrained partitions. 

We turn now to the analysis of unconstrained partitions. The average partition 
function in this case is given by 



(Z U (T)) = J2 ^(T)) 



2T 



dy 



7T l + (2Ty)^ 



*NG u (y) 



where 



G„(y) = ln 



sin (2y) 

2y 



(26) 



(27) 



As before, in the limit of large N the integral over y can be carried out via a 
saddle-point integration. Since the saddle-point is y s = 0, the final result is simply 



(Z U (T))=2 N T X I^- 



(28) 



which yields 



El 



ttN 
6e 2 



0.266 2 



-N 



N. 



(29) 



It is interesting to compare this result with the average energy of a randomly chosen 
configuration s. This quantity, which is defined by 



E r = 2- N il f da, \zZ aiSi I' 

i J ° 8i=±l i 



(30) 



is easily calculated and yields E r = y / 2N/3n for large N. Moreover, comparing 
equations ( p3| ) and (^9|) we find that the lower bound for average optimal energy of 
even partitions (m = 0), which minimizes E^, is larger than that of unconstrained 
partitions by a factor N 1 ^ 2 . The fact that these quantities do not coincide indicates 
that, for unconstrained partitions, m is not a self-averaging quantity, even in the 
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large N limit, i.e. the values of m associated to the best unconstrained partitions 
depend on the specific realization of the set of random variables {%}. In figure 
3(a) we present the numerical estimate for the average optimal energy E u — (E u ) 
obtained through the exhaustive search for N < 24. The data are very well fitted 
by E u = 1.12 2~~ N \/~N . In figure 3(6) we show the ratio r u = ^J~o\jE u as a function 
of N. As before, the finding that this ratio tends to 1 for increasing TV indicates 
that E u is not self-averaging. 



3 Average number of local minima 

As mentioned before, a minimum in the 1-swap neighborhood is a state that has 
lower energy than all the N states that differ from it by a single spin only fllTfl . 
In the statistical mechanics context, these states are usually termed metastable 
states p5[ . The following analysis will be restricted to unconstrained partitions 
only, since for constrained partitions we would have to consider the simultaneous 
flip of two spins in order to satisfy the cardinality constraint. The average number 
local minima with energy E — \ t | is defined by 

(m (t)) =^E 5 f*-E s ^ j n e (i 1 - 2s ^ i-i*D^ ( 3i ) 

where 5 (x) is the Dirac delta function and Q(x) — 1 if x > and otherwise. As 
the calculation is straightforward we will only sketch it in the sequel. Using the 
integral representation of the delta function we obtain 

- e M J] Yl / da > e ~ US3aj (I * - 2s M I - I 1 I) ■ ( 32 ) 

Hence the integral over aj and the summation over Sj can readily be performed, 
yielding 

f°° dt -I e~ %u - e l * -I- e lt - 1 \ 
(M (*))=/ ^e*» UE=\t\ <1, (33) 



and 



2tt \ it 



(■M (*)> = / ^ ^ (^r) = if ^ =1 * I > 1, (34) 



where we have used the interesting result that the integral in equation (|34| ) vanishes 
for all N Pq| . Thus, there are no local minima with E > 1. As usual, for large 
./V the integral in equation (|33| ) can evaluated via a saddle-point integration. The 
final result is 
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where 

and = -d 2 H(£)/d£ 2 . Here, £ is the solution of 

2 2cosh£-e _ ^t 



(36) 



£ sinh^-e-*«/ 2 sinh(^/2) 



= 0. (37) 



The function if (£) is a monotonically decreasing function of E =| i |. In particular, 
it decreases from In 2 at £7 = (£ = 0) to -co at E 1 = 1 (£ = — oo). It vanishes 
at E w 0.851, so the average number of local minima with energy larger than that 
value decreases exponentially with N. 

A more interesting quantity is the average number of local minima regardless 
their energy values, which is defined by 

/oc 
dt(M(t)). (38) 
-oo 

From the above discussion, it is clear that only the close neighborhood of t = 
contributes to this integral, so we can expand the integrand of ( |33|) in powers of t 
and t and keep the lowest order terms only. The final result is 

2.764 w . (39) 

It is interesting to estimate the dependence on N of the typical energy of a local 
minimum. This quantity, denoted by E t , is defined by 

Jdtt M(t) s . , 

Et - <7*W } ' (40) 

which, in the annealed approximation framework p5[ , is approximated by 

Jdtt(M(t)) 
Et (M> ■ (41) 

The procedure to evaluate (fll|) is identical to that used in the evaluation of ( ^8| ) 
and yields 

Ek*—. (42) 
N K ' 

We note that while equation (^) gives the exact leading order term of the average 

number of local minima, equation ( f42| ) is an uncontrolled estimate for the energy of 

a typical minimum. These quantities can be easily estimated numerically: for each 

value of N, ranging from 100 to 3000, we generate 10 5 random states s and count the 

fraction of them that are local minima and measure their energies. We find that the 

numerical data are very well fitted by the equations (M) m (2.81 ± 0.02) 2 N /N 3 / 2 

and E t ~ (1.76 ± 0.04)/iV, which are in quite good agreement with the theoretical 

predictions. 
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4 Conclusion 



To appreciate some of the drastic features of the energy landscape associated to 
the NPP or, equivalently, to the random anti-ferromagnetic Ising model defined by 
the Hamiltonian (^), we compare our results with those of the SK model, which is 
defined by the Hamiltonian 

t j>i 

where the couplings Jy are Gaussian statistically independent random variables 
of zero mean and variance 1/N. In this model the annealed lower bound for the 
ground state energy is E a = — 0.833TV and the number of metastable states 
increases as e 0199N (^8). Hence, in the NPP there are much more local minima 
and the global minima are much deeper than in the SK model. These findings may 
explain the failure of local search techniques to produce good solutions to the NPP. 

Some comments regarding the comparison of our approach with that of Kar- 
makar et al. |l8[ ] are in order. Those authors have derived bounds on the probability 
of occurrence of the event M(E) = 0, where Af(E) stands for the number of states 
s with energy smaller than E. Interestingly, these bounds are related to the first 
two moments of TV: 

(TV 2 ) - (TV) 2 

1 - (TV) < Pr{N = 0} < {J ^ } . (44) 

The first inequality follows trivially from the fact that TV > 0, while the second is an 
improvement of Chebyshev's inequality. Only unconstrained and even partitions 
(m = 0) were considered. However, as acknowledged by Karmakar et al. ]l8| , 
these bounds give no information on the average value of the difference for the best 
partition, except perhaps for its scaling with TV. Also, we should mention that 
Fu has actually carried out a replica analysis of the NPP for the case of even 
partitions to = 0. However, since in that analysis it is assumed that Eq is extensive, 
it misses the low-temperature phase completely. 

Acknowledgments JFF thanks Pablo Moscato for useful conversations. This 
work was supported in part by Conselho Nacional de Desenvolvimento Cientffico c 
Tecnologico (CNPq). 
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Appendix 



In this appendix we calculate exactly the average optimal energy in the regime 
where E m scales linearly with N . Similarly to equation ( |l0|) we write the partition 
function defined in ([|) as 



— oo •/ — oo 



2?r \_ 2tt 



J]] enp[—isj (ajX + rh)] , (45) 

j s j= ±l 

where /3 = 1/T is the inverse temperature. As in the annealed approximation, the 
summation over Sj can easily be performed, yielding 

Z m (T) = N/3 2 dx ^- — exp [-N/3 (xx+ \ x \ +mm)] 



x exp 



- ioo 27Ti J-iir/0 27r * 

(46) 



N da In 2 cosh /3 (aia + m) 
. Jo 

where we have used the self-averaging property 
1 f 1 

— In 2 cosh (3 (xaj + m) = J da In 2 cosh (3 (xa + to) , (47) 



which is exact for N — > oo. In this limit we can carry out the integrals using 
the saddle-point method, and so we obtain the following equation for the average 
free-energy density f m = F m /N: 



1 - 1 



f m = xx+ | x I +mm — — / da In 2 cosh f3 (xa + m) . (48) 
P Jo 

In the zero-temperature limit (/3 — > oo), the saddle-point equations yield x = — 1, 
to = (1 + m)/2 and 

where we have assumed x > and to > 0. The average optimal energy is obtained 
by taking the zero-temperature limit in equation (^) which yields f m — ► E m /N =\ 
x\. 



11 



References 

Anderson P W 1986 Physica A 140 405 

Anderson P W and Fu Y 1986 J. Phys. A: Math. Gen. 19 1605 
Mezard M and Parisi G 1986 J. Phys. France 47 1285 

Cerf N J, Monvel J B, Bohigas O, Martin O C and Pcrcus A G 1997 J. Phys. 
I France 7 117 

Korutcheva E, Opper M and Lopez L 1994 J. Phys. A: Math. Gen. 27 L645 
Fontanari J F 1995 J. Phys. A: Math. Gen. 28 4751 
Inoue J 1997 J. Phys. A: Math. Gen. 30 1047 
Kirkpatrick S and Selman B 1994 Science 264 1297 
Monasson R and Zecchina R 1996 Phys. Rev. Lett. 76 3881 
Hayes B 1997 Am. Scient. 85 108 

Mezard M, Parisi G and Virasoro M A 1987 Spin Glass Theory and Beyond 
(Singapore: World Scientific) 

Garey M R and Johnson D S 1979 Computers and Intractability: A Guide to 
the Theory of NP- Completeness (San Francisco: Freeman) 

Fu Y T 1989 Lectures in the Sciences of Complexity ed D Stein (Reading, MA: 
Addison- Wesley; London: Longman) 

Shamir A 1979 Proceedings 11th Annual ACM Symposium on Theory of Com- 
puting pp. 118 

Kirkpatrick S, Gelatt C D and Vecchi M P 1983 Science 220 671 

Diaz B A and Moscato P 1996 Las Nuevas Tecnicas Heuristicas y las Redes 
Neuronales ed B A Diaz (Madrid: Paraninfo) 

Johnson D S, Aragon C R, McGeoch L A and Schevon C 1991 Operations 
Research 39 378 

Karmakar N, Karp R M, Lucker G S and Odlyzko A M 1986 J. Appl. Prob. 
23 626 

Karmakar N and Karp R M 1982, The Differencing Method of Set Partition- 
ing, Report No. UCB/CSD 81/113, Computer Science Division, University of 
California, Berkeley 

Ruml W, Ngo J T, Marks J and Shieber 1996 Journal Optimization Theory 
and Applications 89, 251 

Toulouse G and Vannimenus J 1980 Phys. Rep. 67 47 



12 



[22] Vannimenus J and Mezard M 1984 J. Phys. Lett. 45 LI 145 

[23] Feller W 1957 An Introduction to Probability Theory and its Applications, vol. 
II (New York: Wiley) 

[24] Daniels H E 1954 Ann. Math. Stat. 25 631 

[25] Tanaka F and Edwards S F 1980 J. Phys. F: Metal Phys. 13 2769 
[26] Erdclyi A 1954 Tables of Integrals Transforms (New York: McGraw-Hill) 
[27] Sherrington D and Kirkpatrick S 1975 Phys. Rev. Lett. 35 1972 
[28] Bray A J and M A Moore 1980 J. Phys. C: Solid St. Phys. 13 L469 



13 



Figure captions 



Fig. 1(a) Average optimal energy obtained through exhaustive search as a function 
of N for even partitions (m = 0). 

Fig. 1(b) Ratio between the standard deviation and the average value of the 
optimal energy obtained through exhaustive search as a function of N for even 
partitions (m = 0). 

Fig. 2 Average optimal energy obtained through exhaustive search as a function 
of m for N = 24. The full curve in the principal graph is the theoretical estimate 
( p5| ) while the one in the inset is the annealed lower bound (pl|). 

Fig. 3(a) Average optimal energy obtained through exhaustive search as a function 
of N 1 / 2 for unconstrained partitions. 

Fig. 3(b) Same as figure 1(6) but for unconstrained partitions. 
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figure 1(a) 
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figure 1(b) 
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figure 3(b) 



